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ABSTRACT 

We use large volume N-body simulations to predict the clustering of dark matter 
in redshift space in f{R) modified gravity cosmologies. This is the first time that the 
nonlinear matter and velocity fields have been resolved to such a high level of accuracy 
over a broad range of scales in this class of models. We find significant deviations from 
the clustering signal in standard gravity, with an enhanced boost in power on large 
scales and stronger damping on small scales in the f(R) models compared to GR at 
redshifts z < 1. We measure the velocity divergence {Pee) and matter (Pss) power 
spectra and find a large deviation in the ratios y/ Pgg / Pss and Pse/Pss between the 
f(R) models and GR for 0.03 < fe/(/i/Mpc) < 0.5. In linear theory these ratios equal 
the growth rate of structure on large scales. Our results show that the simulated ratios 
agree with the growth rate for each cosmology (which is scale dependent in the case 
of modified gravity) only for extremely large scales, k < 0.06/i/Mpc at z = 0. The 
velocity power spectrum is substantially different in the f(R) models compared to GR, 
suggesting that this observable is a sensitive probe of modified gravity. We demonstrate 
how to extract the matter and velocity power spectra from the 2D redshift space power 
spectrum, P(k, fx), and can recover the nonlinear matter power spectrum to within 
a few percent for k < 0.1/i/Mpc. However, the model fails to describe the shape of 
the 2D power spectrum demonstrating that an improved model is necessary in order 
to reconstruct the velocity power spectrum accurately. The same model can match 
the monopole moment to within 3% for GR and 10% for the f(R) cosmology at 
k < 0.2/i/Mpc at z = 1. Our results suggest that the extraction of the velocity power 
spectrum from future galaxy surveys is a promising method to constrain deviations 
from GR. 

Key words: Methods: iV-body simulations - Cosmology: theory - large-scale struc- 
ture of the Universe - dark energy - Modified gravity 



1 INTRODUCTION 

The clustering of galaxies on different scales is a key ob- 
servational tool in the quest to exp lain the current ac- 
celerating expansi o n of the Universe (IPercival et al. 20071: 



Guzzo et al] 120081: iBlake et al\ l20ld. 12011 
20071: ISanchez and PES Collaboration! l20ld; 



2009 



2009 



Schlcgcl et al 



Schlegel et al 



.Laureiis et al] |2011| ; ILSST Science Collaborations 
iGreen et al]\20lA The accelerating expansion may 
be the result of a dark energy component which behaves 
as a repulsive form of gravity or it may be that Einstein's 
theo ry of gravity breaks down on cosmological scales (see 
e.g. iBertschinger fc Zukinll200g| ; IWeinberg et qZ.||2012h . For 



a given cosmology with a smooth dark energy component, a 
measurement of the expansion history gives a prediction for 
the growth rate of structure. Independent measurements of 
the growth rate can be obtained by measuring the clustering 
of galaxies in redshift space, where peculiar velocities distort 
the clustering signal along the line of sight. By testing the 
consistency between the measured growth rate and the pre- 
diction from the expansion history, it is possible to constrain 
models of modified gravity and to disting uish them from a 
smoo t h dark energy compone nt (see e.g. iMortonson et all 
120091 : IVanderveld et al] |2012h . In this paper we measure 
the anisotropic power spectrum in redshift space from large 
volume N-body simulations of f(R) modified gravity and 
general relativity (GR) cosmologies. 
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The f(R) class of models can mimic the effect of a cos- 
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Figure 1. The ratio of the linear growth rate in the F4 (blue), 
F5 (red) and F6 (black) cosmologies to that in ACDM. The ra- 
tios are shown at z = (solid), z = 0.4 (dashed) and z = 1 
(dot dashed). (See Section 2.2 for the description of the modified 
gravity models.) 



mological constant and is set up by modifying the Einstein- 
Hilbert acti on with an arbitrary function of the Ricci sca lar, 
R (see e.g. ICarroIl et aH |200a iNoiri fc Odintsovll2003l ). A 
key feature of these models is the existence of a 'fifth force', 
due to an extra propagating scalar field. Departures from 
general relativity on s mall scales are highly constrained by 
solar system tests (e.g. Will 2006). As a result viable f(R) 
theories must ex hibit a screening mechanis m, the so-called 
chameleon effect jKhourv fc Weltm an 2004), whereby stan- 
dard gravity is recovered in high density environments. The 
range of the fifth force depends nonlinearly on the local cur- 
vature and as a result will change with redshift. The im- 
pact of the chameleon mechanism on the matter and velocity 
fields can only be fully investigated using JV-body simula- 
tions. Any deviations from standard gravity will depend on 
the choice of the function f(R) and the parameter values 
adopted. 

In thi s paper we c o nsider the f(R) model pro- 
posed by iHu fc Sawickil (|2007l ). This modified grav- 
ity model has been incorporated into iV-bod y sim- 
ulations and studied by_ sev eral authors ( Ovaizul 



20081; lOvaizu Lima fc Hul 120081; ISchmidt et all 12 009a : 



Schmidt et alj|2009bl: lLombriser et a/.ll2010l; iFerraro et all 



20ld; Zhao, Li fc Kovamal l201ll ; iLi fc Hul l201ll : 



Gil-Mann et all 



20 111 ) 



A variety of computational box 
sizes from 64 to 400 Mpc/ft on a side have been used. In 
this work we make use of large volume, L hm — 1.5Gpc//i 
and lGpc//i, modified gravi ty simulations using the N- 
body code of ILi et all |2012j ) . The large volume of these 
simulations allow us to study the impact of unique features 
of modified gravity, such as the scale dependent enhanced 
gravitational force, on the clustering signal in redshift 
space. 

Galaxy redshift surveys allow us to study the 3D spa- 
tial distribution of galaxies and clusters. In addition to the 
Hubble flow, galaxies have peculiar velocities, due to lo- 
cal inhomogenities in the density field, which distort the 
measured distances. Measuring the anisotropic distortions 



in the galaxy clustering pattern in redshift space constrains 
/3 = f/b, where 6 is the galaxy bias factor and / is the 
logarithmic derivative of the linear growth rate of struc- 
ture, which is scale independent in the case of gene r al rel - 
ativity. This effect was first described by iKaiserl (|l987l ) 
using linear perturbation theory where the matter power 
spectrum in redshift space can be expressed as a function 
of the power spectrum in real space and /3. Several au- 
thors have extended this linear model to quasi linear scales 
by e . g. including nonlinear velocity terms dScoccimarro 



20041 ; lMatsubarall200Sl ; IPercival fc Whitel2009l ; iTaruva et al 



ZU1U J or Dy considering a phase space distribution tunction 
approach iSeliak fc McDonald! I2011T ) or into the nonlinear 
regime by including the contribut i on of peculiar velocities 
on small scales (e.g. iPeebles |[l976l ; iPeacock fc Doddslll994l ; 
iReid fc Whitell201ll ). 

Modelling the clustering of the dark matter and galax- 
ies in redshift space is extremely challenging. Most models 
contain free parameters such as the linear bias, which quanti- 
fies the difference in clustering between the dark matter and 
galaxies on large scales, and the velocity d ispersion due to in- 
coher ent motions on small scales (see e.g. lOkumura fc Jingl 
l201ll ). These parameters must be included when fitting any 
model and can weaken the constraints on the growth rate. 

Many redshift space distortion models which are cur- 
rently used are only accurate for a limited rang e of scales or 
for g a laxies with a particu la r linear bias ( e.g. | Scoccimarro] 
120041 ; IReid fc White) l201ll ). IScoccimarrd (|2004l ) proposed 
a simple quasi-linear model which includes the nonlin- 
ear velocity power spectrum. B y comparing with measure- 
ments f rom JV-body simulations. I Jennings. Baugh fc Pascolil 
(|201 lal lbh showed that this model performed better than 
commonly used models (which we discuss in Section 3) 
and is accurate on scales k < 0.3/i/Mpc and can re- 
cover the linear growth rate to within a few percent. The 
nonlinear velocity terms in this model may be obtained 
using either a fitting formula calibrated against iV-body 
simulations or from perturbation theory |Matsubaral 1200a ; 
I Jennings. Baugh fc Pascolill2011al ). 

The WiggleZ Dark Energy Survey (|Blake et aZ.I l201lh 
recently measured the growth rate at z = 0.78 to be 
/ = 0.70T0.08 using redshift-space distortions in the galaxy 
power spectrum. Measurements of the linear growth rate 
are degenerate with the bias or clustering amplitude in the 
power spectra and so contraints on the growth rate are often 
quoted as contraints on fas, where as is the rms variance 
in the linear mat ter power spectra smoo thed in spheres of 
radius 8 Mpc/h (|Percival fc White! 120091 ). The 6dF Galaxy 
Survey modelled the 2D galaxy correlation function and 
obtained a low redshift measurement of the growth rate, 
fas = 0.423 ± 0.0 55, at an effective redshift of z=0.067 
jBeutler et aZ.ll2012t) . Recent measurements from the SDSS 
III BOSS survey found dos /dlna = 0.43 ± .069 at an ef- 
fective redshift of z = 0.57 l|Reid et al.ll2012l ). All of these 
results are consistent with the ACDM model and standard 
gravity. Current surveys do not have sufficient precision to 
rule out viable modified gravity models such as the f(R) 
models considered in this paper. Future gala xy redshift sur- 
veys, such as the ESA's EUCLID mission (|Laureiis et all 
1201 ll ) and the gr ound-based stage IV dark energy experi- 
ment, BigBOSS (|Schlegel et al.ll2009h . aim to measure the 
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growth rate to within 2%, which will place significant con- 
straints on currently allowed modified gravity models. 

A key feature of redshift space distortion models is that 
the linear growth rate is assumed to be scale independent. 
This assumption is not true for the modified gravity cosmol- 
ogy considered in this paper (see Fig.Q]). In addition, several 
models of redshift space distortions suffer from systematic 
biases when fitting for a scale in dependent growth rate over 
a rang e of scales (see Figure 5 in Ijennings. Baugh fc Pascolil 
l2011bf ). In order to avoid assuming a specific scale depen- 
dence for the growth rate, we instead focus on recovering the 
velocity and matter power spectra as a function of scale us- 
ing the full two dimensional redshift space power spectrum. 
This approach makes use of the full 2D power spectrum mea- 
sured from a survey and the extracted matter and velocity 
power spectra could be compared to predictions from the 
standard cosmological model. 

In this paper we measure the power spectrum in red- 
shift space fr om a suite of l arge volume simulations of f(R) 
cosmologies |Li et a/.ll20"l2l ). This is the first time that pre- 
dictions for the redshift space clustering in f(R) modified 
gravity models have been presented. The resolution of our 
simulations allows us to accurately resolve the nonlinear 
matter and velocity fields and quantify the deviations from 
a model of general relativity. Here, we restrict our study to 
the clustering of the dark matter. We examine the differ- 
ence between the velocity power spectra in each cosmology 
and its importance in modelling the redshift space cluster- 
ing signal in both the standard and modified gravity model. 
We test how well quasi linear models for the redshift space 
distortions describe the amplitude and shape of the mea- 
sured power spectrum. A follow up paper will examine the 
redshift space distortions in the clustering of halos as well as 
testing n onlinear models for redshift space distortions (see 
also e.g. IMaruUi et oEll2012i . for a recent study of redshift 
space distortions in interacting dark energy models). 

This paper is organised as follows: In Section [2] we dis- 
cuss the f(R) modified gravity cosmological model and de- 
scribe the TV-body simulations used in this paper. In Sec- 
tion0we review the theory of redshift space distortions and 
present the models which will be tested. The main results 
of the paper are presented in Section [4] Measurements of 
the redshift space power spectra for both general relativ- 
ity and the f(R) models are presented in Section 14.11 In 
Section 14.21 we present the velocity power spectrum mea- 
sured from the simulations. Using a quasi linear model for 
the redshift space power spectrum we attempt to extract 
both the matter and velocity power spectra from the two 
dimensional redshift space power spectrum in Section 14.31 
In Section 14.41 we examine how well the moments of the 
redshift space power spectrum can be recovered using this 
quasi linear model. Our conclusions are presented in Sec- 
tion [5] Throughout the paper we shall use the unit c = 1 
and metric convention (+,—,—,—). Greek letters p,u,--- 
run over 0, 1, 2, 3 and Latin letters i, j, k, ■ ■ ■ run over 1,2,3. 



2 F{R) COSMOLOGIES 

This section gives the theoretical background for the modi- 
fied gravity model considered in this paper. We outline f(R) 
cosmologies in Section 2.1, explain the chameleon mecha- 



nism in Section 2.2 and describe the JV-body code and sim- 
ulations in Section 2.3. 



2.1 The f(R) gravity model 

The f{R) gravity model is a straightforward generalisation 
of GR: the Ricci scalar, R, in the Einstein-Hilbert action, 
S, is replaced with an algebraic function, f(R) (see e.g., 
ISotiriou fc Faraonil l20ld : |Pe Felice fc Tsuiikawal |2010l . for 
recent reviews): 

S = j A A xsr-gl^[R + f{R)]+C m \, (1) 

in which Mpi is the Planck mass, Mpj 2 = 8nG with G being 
Newton's constant, g is the determinant of the metric <; M „ 
and C m is the Lagrangian density for matter fields (photons, 
neutrinos, baryons and cold dark matter). By specifying the 
functional form of f(R) one specifies the ,f(R) gravity model. 

Varying the action defined in Eq. ([1} with respect to 
the metric p M „ yields the modified Einstein equation 

G M „ + f R R^ - Q/ - D/b) g^ - V M V„/ fl = 8ttGT™ , (2) 

in which G M ,y = R^v — ^g^ v R is the Einstein tensor, 
f R = df/dR, V M the covariant derivative compatible with 
the metric <? M „, □ = V a V a and T^J is the energy momentum 
tensor for matter. One can consider Eq. ((2]) as a fourth-order 
differential equation, or alternatively the standard second- 
order equation of GR with a new dynamical degree of free- 
dom, /ij, the equation of motion of which can be obtained 
by taking the trace of Eq. @ 

af R =^(R-f R R + 2f + 8TvGp m ), (3) 

where p m is the matter density. The new degree of freedom 
fn is sometimes dubbed the scalaron in the literature. 

Assuming that the background Universe is described 
by the flat Friedmann- Robertson- Walker (FRW) metric, the 
line element in the perturbed Universe is written as 

ds 2 = a{rf) [(1 + 2$)d77 2 - (1 - 2*)da; i dx i ] , (4) 

in which r\ and x % are, respectively, the conformal time and 
comoving coordinates, 4>(ry, x) and ^(77, x) are the Newto- 
nian potential and perturbation to the spatial curvature, 
and are functions of both time (r\) and space (x); a denotes 
the scale factor of the Universe where a — 1 today. 

As we are mainly interested in the large-scale structures 
much smaller than the Hubble scale, and since the time vari- 
ation of fa is very small in the models considered below, we 
shall work in the quasi-static limit by neglecting the time 
derivatives of f R . In this limit, the scalaron equation, Eq.[2] 
reduces to 

V 2 / fl = -ia 2 [R{f a )- R+8irG(p m - p m )] , (5) 

in which V is the three dimensional gradient operator (to be 
distinguished from the V introduced above), and the overbar 
represents the background value of a quantity. Note that R 
can be expressed as a function of f R . 

Similarly, the Poisson equation which governs the New- 
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tonian potential, $, can be simplified to 



V $ = — — a (p„ 



+ -a 2 [R(f R ) 



(6) 



by neglecting terms involving time derivatives, and using 
Eq. ([5| to eliminate V 2 /h. 

According to the above equations, there are two poten- 
tial effects of the scalaron on cosmology: (i) the background 
expansion of the Universe may be modified by the new terms 
in Eq. ((2J and (ii) the relationship between gravity and the 
matter density field is modified, which can change the mat- 
ter clustering and growth of density perturbations. Clearly, 
when \f R \ -C 1, we have R ss —8nGp m from Eq. ([5]) and so 
Eq. © reduces to the normal Poisson equation in GR; when 

is large, we instead have \R — R\ <C 87rG|p m — p m \ and 
so Eq. (fSj reduces to the normal Poisson equation with G 
rescaled by 4/3. Note that this factor of 4/3 is the maximum 
enhancement of gravity in f(R) models, independent of the 
specific functional form of f{R)- The choice of f{R), how- 
ever, is important because it governs when and on which 
scale the enhancement factor changes from unity to 4/3: 
scales much larger than the range of the modification to 
Newtonian gravity mediated by the scalaron are unaffected 
and gravity is not enhanced there, while on much smaller 
scales the 4/3 enhancement is fully realised - this results 
in a scale-dependent modification of gravity and therefore a 
scale-dependent growth rate of structure (see Fig. 1). 

The relationship between $ and is also changed in 
f(R) models, with the remaining components of the modi- 
fied Einstein equation giving 



V 2 (* - $) 



(7) 



where we have assumed that \f R \ -C 1. This implies that 

V 2 ($ + *) = SnG ( Pm - pm) a 2 . (8) 

Therefore the relationship between the lensing potential and 
the matter density perturbations remains unchanged in f(R) 
gravity models. 



2.2 The chameleon mechanism 

The f(R) gravity would be ruled out by local tests of grav- 
ity due to the factor of 4/3 enhancement to the strength 
of Newtonian gravity. Fortunately, it is well known that, 
if f( R) is chosen appropriatel y , the chameleon m echa- 
nism (|Khourv fc Weltmanl 12004 iMota fc Shawl 120071 ) can 
be exploited to suppress the enhancement allowing this 
class of models to satisfy experimental constraints in 
high matter density regions such a s in our Solar sys- 



tem ([Navarro fc Van Acoleven 



iHu fc Sawickill2007l : iBrax et alil2008f ) 



20071: iLi fc Barrowl 120071: 



The essence of the chameleon mechanism is as follows. 
The modifications to the Newtonian gravity can be consid- 
ered as an extra, or fifth force mediated by the scalaron. 
Because the scalaron itself is massive, the force is of the 
Yukawa type and is suppressed by the exponential factor 
~ exp(— mr) in which m is the scalaron mass and r the 
distance between two test masses. In high matter density 
environments, m is very heavy and the suppression becomes 
very strong. In practice, this is equivalent to the fact that 
I /ill -C 1 in high density regions because of the exponential 
suppression, which leads to the GR limit as discussed above. 



As a result, the functional form of f(R) is crucial to 
determine whether the fifth force is sufficiently suppressed 
in high density envi ronments. In th i s wor k we study the f(R) 
model proposed bv lHu fc Sawickil |2007h . for which 



f(R) 



-M' 



Cl 



-r/m 2 y 



c 2 (-R/M 2 ) n + 1 ' 



(9) 



with M 2 = 8nGpmo/3 = H^Qm, where H is the Hubble 
expansion rate and Q m is the present-day fractional density 
of matter. Hereafte r a subscript alway s means the current 
value of a quantity. IHu fc Sawic kl (|2007l) demonstrated that 
I /ho I < 0.1 is required for this model to evade Solar system 
constraints, although the exact value also depends on the 
behaviour of f R in the Galaxy. 

In the background cosmology, the scalaron f R always 
sits at the minimum of the effective potential which governs 
its dynamics, defined as 



Kff CM = - (R - f R R + 2R + 87rGp„ 



(10) 



around which it oscillates quickly (|Brax et al.ll2012h . There- 
fore we have 



R « S-nGpr, 



2/ = 3Af 2 (a" 3 + ^i 

3C2 



(11) 



To match the ACDM background evolution, we need to have 

(12) 



where Qa is the current fractional energy density of the dark 
energy (cosmological constant). 

By taking S1a ~ 0.76 and f2 m ~ 0.24, we find that 
\R\ « 41M 2 3> M 2 , and this simplifies the expression of the 
scalaron to 



Cl 



f R w -n-j — 



M 2 
~R 



n + l 



(13) 



Therefore, two free parameters, n and ci/c 2 ., completely 
specify the f(R) model. Indeed, the latter is related to the 
value of the scalaron today, fuo, as 



Cl 



3(1 + 4 



/ho- 



(14) 



In what follows we shall study three f(R) models with n = 
1 and | /ho | = 10~ 6 , 10~ 5 , 10~ 4 , which will be referred to 
as F6, F5 and F4 hereafter, respectively. These particular 
parameter choices arise from cluster abundance constraints 
on the f(R) gravity m odel. The current constaint found by 
ISchmidt etaU (|2009bl ) is \f R0 \ < l.3±o 7 6 x 10" 4 taking into 
account mass calibration errors. 



2.3 TV-body simulations of f(R) gravity 

From Eqs. ([S] [5]) we can see that, given the matter density 
field, we can solve for the scalaron field, f R , from Eq. ([5]) 
and plug this into the modified Poisson equation ([6]) to solve 
for $. Once $ is at hand, we can difference it to calculate 
the (modified) gravitational force which determines how the 
particles move subsequently. This is exactly what we need 
to do in iV-body simulations to follow the evolution of the 
matter distribution. 

The main challenge in the iV-body simulation of models 
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Figure 2. Top panels: The two dimensional power spectrum measured from the GR simulation as a function of wavenumber perpen- 
dicular, fc pcrp , and parallel, fc para , to the line of sight. The colored shading and lines represent the amplitude of the power spectrum, 
log 10 P as indicated by the labels and the scale bar at the top. The real space power spectrum has been plotted at fc porp — > — fc pcrp to 
allow comparison with the redshift space P(k). Bottom panels: The two dimensional power spectrum measured from the F4 simulation 
as a function of modes perpendicular, fc pclp , and parallel, fc pal - a , to the line of sight. 



such as f(R) gravity is to solve the scalaron equation, Eq. 
which is in general very nonlinear. For this we need a mesh 
(or a set of meshes) on which fn can be solved. This implies 
that mesh-based iV-body codes are the most suitable for 
this task. On the other hand, particle-based codes are more 
difficult to apply in this case, as we do not have an analytical 
formula for the modified force law (i.e. the equivalent of r~ 2 
in the Newtonian case). 

iV-body simulations 
lated theories have been 



for f(R) gr avity 
er formed by Ovaiz 
Q£D; ISchmidt et al 



Ovaizu. Lima fc Hu 

Zhao, Li fc Kovamal d201lf); ILi fc Zhaol J2009I 
Schmidt I (120091); iLi fc Barrow! j20 111 ): iBrax etaU 
Davis et all (|2012h . However, these simulations were af- 




fected by the small box size used and the limited resolution. 
For this work we have run simulations of f (R) cosmologie s 
using the recently developed EC0SMDG code (|Li et aZ.|[2012f l. 
EC0SMQG is a modification of the mesh-based iV-body code 



RAMSES l |Tevssierll2003 ). which calculates the gravitational 
force by first solving the Poisson equation on meshes using 
a relaxation method to get the Newtonian potential and 
then differencing the potential; it does not solve gravity by 
summing over the forces from nearby part icles as in don e 
for example in the simulation code GADGET (|Springej|2005h . 
Additional features of the EC0SM0G code include: 

(i) The adaptive mesh refinement (AMR) , which refines a 
mesh cell, i.e. splits it into 8 son cells, if the number of par- 
ticles in a cell exceeds a pre-defined number (the refinement 
criterion). As such it gives higher force resolution in high 
matter density regions where the chameleon effect is strong 
and the f(R) equation is more nonlinear. The refinement 
criterion is normally chosen as a particle number between 
8 and 12, and in our simulations we adopt a condition of 9 
particles. This adaptive mesh allows us to reach comparable 
spatial resolution to codes like GADGET. 

(ii) The multigrid relaxation algorithm that ensures quick 
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Figure 3. Top panels: The two dimensional power spectrum measured from the F5 simulation. Bottom panels: The two dimensional 
power spectrum measured from the F6 simulation. The coloured shading correspond to log 10 P as shown in the color bar in Fig. [2] As in 
Fig. [2] the left panels show the power spectrum in real space and the right panels show rcdshift space. 



convergence. The relaxation method finds the solution to an 
elliptical partial differential equation (PDE) on a mesh by 
iteratively updating the initial guess until it converges, i.e., 
becomes enough close to the true solution. But the rate of 
convergence slows down quickly after the first few iterations. 
To improve on this, one can 'coarsify'the PDE, i.e., move it 
to a coarser mesh, solve it there and use the coarse solu- 
tion to improve the solution on the original fine mesh. Un- 
like other codes, EC0SM0G does this on all the AMR meshes, 
greatly improving the convergence properties of the whole 
code. 

(iii) The massive parallelisation which makes the compu- 
tation very efficient. This is the key feature that enables us 
to run large simulations such as the ones to be described 
below, which are bey ond the rea c h of any se ri al code, such 
as tho se developed bv lLi fc Zhaol l |2009ll2010h ; lLi fc Barrowl 
|201ll ). 

A convergence criterion is used to determine when the re- 
laxation method has converged. In EC0SMQG convergence is 
considered to be achieved when the residual of the PDE, i.e., 



the difference between the two sides of the PDE, is smaller 
than a predefined parameter e. We have checked that for 
e < 10~ 8 the solution to the PDE no longer changes sig- 
nificantly when e is reduced further. Our choices of e are 
listed in Table 1. The computational time depends on both 
the value of e and the model. The f(R) gravity simulations 
can take a few times longer t o run than the GR simulation. 
More details can be found in iLi et al\ (|2012l ). 



For the study of redshift space distortions, large sim- 
ulations boxes are essential to accurately model behaviour 
on very large scales. For this reason, we have run two sets 
of simulations, with Lbox = 1.0ft _1 Gpc and 1.5/i _1 Gpc re- 
spectively. The initial cond itions are generated at z — 49 
using the MPgrafic code (|MPgrafich . and each suite of 
F4/F5/F6/GR simulations uses the same initial conditions 
because at z — 49 the difference in the matter distribution 
in the different cosmologies is negligible. The specifications 
of the simulations are summarised in Table [1] 
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Table 1. Some technical details of the simulations performed in this work. F6, F5 and F4 are respectively the abbreviations which 
denote the f(R) models with |/ho| = 10 -6 , 10 — 5 , 10 -4 . For all model s we have assumed O m = 0.24 and Q\ = 0.76, and to generate 
the initial conditions we have used <r$ = 0.769 (in agreement with e.g. ISanchez et al. ||2009|) . We use the same initial conditions for all 
models in each simulation, because at the initial time, Zi = 49, t he difference i n the power spectra of different models is negligible, e is 
the residual for the Gauss-Scidcl relaxation used in the code (see lLi et alll2012ri , and the two values of the convergence criterion are for 
the coarsest level and finest levels respectively. Other cosmological parameters are a Hubble constant of Hq = 73km/s/Mpc and a scalar 
spectral index of n = 0.961. 



models ^box particles domain meshes finest meshes convergence criterion realisations 

ACDM, F6, F5, F4 1.0/i _1 Gpc 1024 3 1024 3 65536 s |e| < 10 _12 /10 -8 1 

ACDM, F6, F5, F4 1.5h -1 Gpc 1024 3 1024 3 65536 3 |e| < 10" 12 /10" 8 6 



3 MODELS FOR REDSHIFT SPACE 
DISTORTIONS 

In this section we first review the linear perturbation theory 
of redshift space distortions (Section 3.1) before outlining 
some extended models which go beyond linear theory (Sec- 
tion 3.2). 



3.1 Linear theory 

Inhomogeneous structure in the universe induces peculiar 
motions which distort the clustering pattern measured in 
redshift space on all scales. This effect must be taken into 
account when analysing three dimensional datasets which 
use redshift to estimate the radial coordinate. Redshift space 
effects alter the appearance of the clustering of matter, and 
together with nonlinear evolution and bias, lead the mea- 
sured power spectrum to depart from the simple predictions 
of linear perturbation theory. The comoving distance to a 
galaxy, s, differs from its true distance, x, due to its pecu- 
liar velocity, v(x) (i.e. an additional velocity to the Hubble 
flow). The mapping from redshift space to real space is given 
by 



x + u x z, 



(15) 



where u z = v ■ z/(aH) and H(a) is the Hubble param- 
eter. This assumes that the distortions take place along 
the line of sight denoted by z (N.B. this is the plane 
parallel approximation). This assumption will break down 
for some pai rs of galaxies in a survey which has a wide 
field of view (|Raccan clli 2010). Nevertheless, the impact of 
this systematic on clustering statistics has been shown to 
be small in comparsio n to the effects of nonlinear growth 
(|Samushia et al\\2Q12\ ). 

On small scales, randomised velocities associated with 
the motion of galaxies inside virialised structures reduce the 
power. The dense central regions of galaxy clusters appear 
elongated along the line of sight in redshift space, which 
produ ces the 'fingers of God' effect in redshift survey cone 
plots (| Jackson! 1 1972h . On large scales, coheren t bulk flows 
distort clustering statistics (see lHamilton|[l99c1 . for a review 
of redshift space distortions.). For growing perturbations on 
large scales, the overall effect of redshift space distortions 
is to enhance the clustering amplitude. Any difference in 
the velocity field due to mass flowing from underdense re- 
gions to high density regions will alter the volume element, 
causing an enhancement of the apparent density contrast in 
redshift space, S s (k), compared to that in real space, S r (k). 



This effect was first analyzed bv lKaiserl (|l987h in linear per- 
turbation theory and can be approximated by 



5 B {k) = S r (k)(l + fi 2 /3), 



(16) 



where fj, is the cosine of the angle between the wavevector, 
k, and the line of sight, /3 = f/b, f is the linear growth rate 
and the bias b — 1 for dark matter. 

The 'Kaiser formula' (Eg. I16[l relates the overdensity in 
redshift space to the corresponding value in real space and 
is the result of several approximations, e.g., that the veloc- 
ity and density perturbations satisfy the linear continuity 
equation, 



5 = -f0, 



(17) 
All 



where 9 — V • u is the velocity divergence, 
of these assumptions are valid on scales that are 
well described by linear perturbation theory and will 
break down on differ ent scales as the density fluctu- 
ations grow (see e.g. Ijennings. Baugh fc Pascolil 2011al . 



Scoccimarrd l|2004 ) and 



for more details). As shown in 

Ijennings. Baugh fc Pascolil (^Olla ). the linear regime cor- 
responds to a different range of scales for the matter and 
velocity fields. In particular, linear theory is only a good 
description of the velocity power spectrum on surprisingly 
large scales. We will discuss this further in Section \4. 2 1 

Rather than use the full 2D power spectrum, P(k,fj,), 
it is common to decompose the matter power spectrum in 
redshift space into multi pole moments u sing Legendre poly- 



redshut space into multipole moments ui 
nomials, Li(fi), (see e.g. lHamilton|[l99c^ ) 



(18) 



where the summation is over the order, I, of the multipole. 
The anisotropy in P(k) is symmetric in /i, as P(k, fi) = 
P(k, —fi), so only even values of I are summed over. Each 
multipole moment is given by 
fi 



Pi(k) = 



21 + 1 



P(k,fi)Li(n)dfj,, 



(19) 



where the first two non-zero moments have Legendre poly- 
nomials, Lo (a 4 ) — 1 an d L2(n) = (3/i 2 — l)/2. Using the 
linear model in Eq. 1161 the first three multipole moments 
are given by 

/ l + f/3 + 1/3 2 ^ 
= Pss(k)[ |/3 + |/? 2 




35^ 



(20) 



where Pss(k) = (|<5(fc)| ) denotes the real space matter 
power spectrum. Note we have neglected the superscript s 
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here for clarity. In practice, Pss(k) cannot be obtained di- 
rectly for a real survey wit hout making approximations (e.g. 
iBaueh fc Efstathiodfl994h . 

In this paper we conside r the e stimator for /3 suggested 
by ICole, Fisher fc Weinberg! (|l994l ). which is the ratio of 
the quadrupole to monopole moments of the redshift space 
power spectrum, P2(k)/Po(k), which is independent of the 
real space power spectrum. 

3.2 Quasi-linear and Non-linear models for the 
redshift space power spectrum 

Assuming the line of sight component of the peculiar velocity 
is along the z-axis, the full nonlinear relation between the 
real and redshift space power spectrum c an be written as 
(|Scoccimarro. Couchman fc Friemanlll999f ) 

P S (k,») = J ^3e- lk -(e lAA -[5(x)-e(x)] 

X[«5(x')-0(x')]>, (21) 

where A = fc/x, u z is the comoving peculiar velocity along the 
line of sight, Au z = u z (x) — u z (x'), r = x — x', 9 — V z • u z , 
and the only approximation made is the plane parallel ap- 
proximation. At small scales (as k increases) the exponential 
component damps the power, representing the impact of ran- 
domised velocities inside gravitationally bound structures. 

Simplified models for redshift space distortions are fre- 
quently used. Examples include multiplying Eq. 14.41 by a 
factor which attempts to take into account small scale ef- 
fects, invoking either a Gaussian or exponentia l distribution 
of peculiar velocities l|Peacock fc Doddsl I1994T ). A popular 
phenomenological example of this which incorporates the 
damping effect of velocity d ispersion on small scales is the 
so-called 'dispersion model' (jPeacock fc Doddslll994l l. 

F(M) = ^)(l + ^f lrT ^, (22) 

where a p is the pairwise velocity dispersion along the line 
of sight, which is treated as a parameter to be fitted to the 
data. 

The linear model for the redshift space power spectrum 
can be extended by keeping the nonlinear velocity power 
spectra terms in Eq. 1211 The velocity divergence auto power 
spectrum is the ensemble average, Peg = (\9\ 2 ) where 9 = 
V • u is the velocity divergence. The cross power spectrum of 
the velocity diverg ence and matter density is Pse = (\S9\). 
IScoccimarrJ (|2004n proposed the following model for the 
redshift space power spectrum in terms of Pss , the nonlinear 
matter power spectrum, Pee and Pse, 

P s (k,v)= (23) 
(Pss(k) + 2[i 2 P se (k) + ^ 4 P e e(k)) x e -( fe ^> 2 , 
where er„ is the ID linear velocity dispersion given by 

2 1 f Pee(k) 3 , . 

ff « = 3 J ~l^ d ^ {M) 

In linear theory, Pee and Pse take the same form 
as Pss and depart from this at different scales. Us- 
ing a simul ation with 512 3 particles in a b ox of length 
479/i~ 1 Mpc (|Yoshida. Sheth fc Diafericll200ll ). IScoccimarrol 
|2004l ) showed that this simple ansatz for P 3 (k, /i) was an 



improvement over the Kaiser formula when comparing to 
the results of TV-body simulations in a ACDM cosmology. 

In nonlinear models for the power spectrum in redshift 
space there is a degeneracy between the nonlinear bias, the 
difference between the clustering of dark matter and halos 
or galaxies, and the scale dependent damping due to velocity 
distortions on small scales. This degeneracy will complicate 
any measurement of the growth rate using redshift space 
clustering information on small scales. In addition these 
models assume that the growth rate is scale independent 
as is the case in general relativity. As shown in Fig. [TJ the 
f(R) model has scale dependent growth rates which would 
have to be included when fitting any model to the measured 
redshift space power spectrum. 

In this paper we analyze the redshift space clustering of 
the dark matter in f(R) and ACDM cosmologies. We restrict 
our analysis to the linear and quasi-linear regime where the 
bias is typically assumed to be scale independent and so our 
results can be more easily extended to linearly biased t rac- 
crs of the dark matter field (but see lAngulo et a/1l2008l . for 
counterexamples). We will model the redshift space distor- 
tions in the clustering of halos in both of these cosmologies 
in future work. 

We shall use the following model for the 2D redshift 
space power spectrum, 

P(k,n) = P ss (k) + 2[i 2 Pse(k)+ ^ 4 Pee{k), (25) 
where the first two multipole moments are given by 

(SS)-(ill)O)- « 

This model has been shown to be a good fit to the power 
spectr um in redshift space measured fr om simulations at 
z < 1 (|jennings. Baugh fc Pascolll2011al rbl). 



4 RESULTS 

In Section l4~P] we present measurements of the redshift space 
power spectra for both general relativity and the f(R) mod- 
els. In Section 14.21 we present the velocity power spectrum 
measured from the simulations. We attempt to extract both 
the matter and velocity power spectra from the two dimen- 
sional redshift space power spectrum using a quasi-linear 
model for the redshift space P(k) in Section 14.31 In Sec- 
tion 14.41 we examine how well the moments of the redshift 
space power spectrum can be recovered using this quasi lin- 
ear model. 

4.1 The power spectrum in redshift space 

In Figs. [2] and [3] we plot the two dimensional power spec- 
trum measured at z = from the GR and F4 simulations 
and the F5 and F6 simulations respectively as a function of 
wavenumber perpendicular, fc porp , and parallel, fc para , to the 
line of sight. The color contours and lines represent the am- 
plitude of the power spectrum, log 10 P. In each figure the real 
space power spectrum has been plotted as fc porp — > — k pe r P 
in order to allow a side-by-side comparison with the redshift 
space P(k). These figures clearly show that the spherical 
symmetry seen in the real space power spectrum (left panel) 
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Figure 4. Left panel: The ratio of the monopole of the redshift space power spectrum to the real space power spectrum at z = in 
the F4 (blue) and GR (black) cosmologies. We plot the ratio measured from two simulation boxes: Lbox = 1000 Mpc/h (circles) and 
Lbox = 1500 Mpc//i (dashed lines). The linear theory prediction for each model is shown as a green dot dashed line for the f(R) model 
and a red dashed line for GR. The shaded regions represent the errors on the ratios measured from six realisations of the f(R) (blue 
hatched) and GR (grey solid) cosmologies. Right panel: The ratio of the quadrupole to monopole moment of the redshift space P(k) at 
z = in the F4 (blue) and GR (black) cosmologies. 
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Figure 5. Left panel: The ratio of the quadrupole to monopole moment of the redshift space power spectrum at z = in the F5 (blue) 
and GR (black) cosmologies. Right panel: The same ratio of the quadrupole to monopole moment of the redshift space power spectrum 
at z = in the F6 (blue) and GR (black) cosmologies. The linear theory prediction for each model is shown as a green dot dashed line 
for the f(R) model and a red dashed line for GR. The shaded regions represent the errors on the ratios measured from six simulations 
of the f(R) (blue hatched) and GR (grey solid) cosmologies. The simulation box size used was Lbox = 1500Mpc//i on a side. 



is distorted in redshift space (right panel): on large scales 
(k — > 0), the amplitude of the redshift space power spectrum 
is increased compared to that in real space, whereas on small 
scales, the power spectrum is damped and elongated along 
the line of sight in redshift space compared to real space. 
These effects were first convincingly observ ed in the clus- 
terin g of galaxies measured by the 2dFGRS (|Peacock et al. I 
l200ll ) and again recently by the SDSS-III BOSS survey and 
the W iggleZ Dark Energy Survey, c.f . Fig. 3 in iReid et aU 
(|2012l ) and Fig. 2 in lBlake et ail (120111 ). These two effects are 
more pronounced in the F4 simulation where the large scale 
boost and small scale damping appear larger than in GR. 
Overall, the redshift space P(k) for the f(R) model appears 



far more distorted and asymmetrical than the corresponding 
redshift space P(k) in GR. 

In the left panel of Fig. [4] we plot the ratio of the 
monopole of the redshift space power spectrum to the real 
space power spectrum at z = measured from the GR 
(black) and F4 (blue) simulations. The right panel of this 
figure shows the ratio of the quadrupole to monopole mo- 
ment of the redshift space power spectrum for both models 
at z — 0. The redshift space power spectra are obtained from 
the simulations after averaging over the P(k) obtained by 
treating the x, y and z directions in turn as the lines of sight. 
The errors on the ratios are plotted as a blue hatched region 
for the f{R) cosmology and as a solid grey region for GR 
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and represent the scatter amongst six realisations. Note we 
have compared the errors obtai ned from six simulations to 
those fr om ten simulations from I Jennings. Baugh fc Pascolil 
|2011bl ). which have the same box size and particle number, 
and find a 20% decrease in the z = error on the largest 
scale, fc ~ 0.01/i/Mpc, when we use six simulations instead 
of ten. The linear theory predictions, given by Eq. 14.41 are 
shown as a green dot dashed line for the f{R) model and 
a red dashed line for GR. These predictions use the linear 
growth rate for each model which for GR is f(z = 0) = 0.42 
and is a scale dependent factor, f(k), for the f(R) cosmology 
(see Fig. 1). 

It is clear from Fig.Uthat the redshift space power spec- 
trum in the f{R) model has a different shape compared to 
that in GR. Firstly, there is an increased boost in the cluster- 
ing signal on large scales, k < 0.07/i/Mpc, due to increased 
bulk flows into overdense regions seen in the modified gravity 
simulation. Nevertheless, the measurement from the simu- 
lation is significantly below the linear theory prediction for 
this cosmology. The small scale damping due to incoher- 
ent motions within virialised structures is also more pro- 
nounced in the f(R) cosmology compared to GR on scales 
k > 0.1/i/Mpc. It is clear from this plot that the linear per- 
turbation theory limit is only attained on extremely large 
scales (fc < 0.02/i/Mpc _ ) for each model. On these scales 
the two models cannot be distinguished within the error 
bars. The large scale boost in the redshift space power spec- 
trum in the f(R) cosmology compared to GR is less pro- 
nounced in the quadrupole to monopole moment ratio plot- 
ted in the right panel of Fig. 3] The increased damping of the 
redshift space power spectrum found in the modified gravity 
model can again be seen on small scales (fc > 0.1/i/Mpc). 
The increased damping of the redshift space power spectrum 
on small scales is a generic feature of f(R) gravity. This is 
because the fifth force causes the particles to move faster, 
i.e., with a larger velocity dispersion than they do in GR for 
the same mass structure. 

In Fig.|4]we also compare the measured ratios from two 
different simulation boxes of L box — 1000 Mpc//i (circles) 
and ibox = 1500 Mpc/h (dashed lines) on a side. Each sim- 
ulation has the same number of particles, 1024 3 , but with 
different resolutions. In the L b0 x = 1500 Mpc/h simulation, 
the domain gricQ has 1024 3 cells with a refinement criterion 
equal to 9 particles (which means that a cell is refined into 
eight 'son' cells if the number of particles inside it exceeds 
9). The Lbox = 1000 Mpc//i simulation was run with the 
same domain grid and refinement criterion. As a result, both 
the mass and force resolutions are higher in the _L b0 x = 1000 
Mpc/h simulation. The force resolution is 15.3 kpc/h and 
22.9 kpc/h in the 1 Gpc/h and 1.5 Gpc/h box simulations 
respectively. This can be compared to the force resolution 
of 781 kpc/h for the 4 00 Mpc/fe box and 125 kpc/h for the 
64 Mpc/h box used bv lSchmidt et al\ i|2009bh . In Fig. H the 
higher resolution simulation shows more damping on small 
scales, and this is clearly seen in the quadrupole to monopole 
moment plotted in the right panel. This difference in the re- 
sults from different computational boxes shows that large 
volume high resolution simulations are essential in order to 



1 Here the domain grid is the finest grid that is uniform across 
the computational box. 
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Figure 6. The ratio of the quadrupole to monopole moment of 
the redshift space power spectrum for ACDM and the F4 model 
at z = 0.4 (bottom panel) and 2 = 1 (top panel). The linear 
theory prediction at each redshift is shown as a green dot dashed 
line for the f(R) model and a red dashed line for GR. 



accurately resolve the velocity field on scales fc > 0.3/i/Mpc 
and provide accurate predictions for the power spectrum 
in redshift space. In this paper we restrict our study to 
scales fc < 0.3ft/Mpc where the velocity field is accurately 
resolved in both simulations. We will only show results from 
the Lbox = 1500 Mpc/h simulation in the rest of the paper 
unless otherwise stated. 

In Fig. [S] we plot the ratio of the quadrupole to 
monopole moment of the redshift space power spectrum at 
2 = measured in the F5 (left panel) and F6 simulations 
(right panel) along with the measured ratios in GR. The fifth 
force is more strongly suppressed in both of these models 
compared to the F4 model. As a result we find no detectable 
boost in the power spectrum on large scales compared to GR 
and less damping on small scales compared to F4. The F5 
model in the left panel of Fig.[5]still shows significantly more 
damping on small scales compared to GR, whereas the dif- 
ference between the F6 model and GR at z = is very small 
even down to small scales (fc ~ 0.3/i/Mpc). This implies that 
the particles inside halos are not significantly affected by the 
fifth force in the F6 model, while in F5 they have started to 
feel an effect. 
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Fig. [6] shows how the quadrupole to monopole moment 
ratios of the redshift space power spectrum for GR and the 
F4 model change with redshift. The lower (upper) panel 
shows the measured ratios at z = 0.4 (z = f.0), together 
with the linear theory predictions for each model (red dotted 
line for GR and green dot-dashed line for F4) at the same 
redshift. We can clearly see that the linear theory predictions 
agree with the measured ratios on slightly smaller scales 
(down to k ~ 0.05/i/Mpc) at z = 1 than z = for both 
models, as expected. The ratio P2/P0 in the F4 models is 
slightly larger than GR on large scales but it suffers stronger 
damping on small scales compared to GR, such that the ratio 
becomes smaller than GR on nonlinear scales. This large 
scale boost of the power spectrum and extra small scale 
damping compared to GR are more pronounced at higher 
redshifts, which seems to contradict the naive expectation 
given that the fifth force is weaker then. This is because 
the fifth force effect has already been felt well before 2 = 1 
as the screening is weaker at these epochs, after which the 
GR result slightly catches up and so the difference from F4 
is reduced. Indeed, the same trend can also be seen in the 
linear perturbation results plotted in this figure. 

Fig. [7]shows the ratio of P2 / Po in the F4 models to that 
in GR at 2 = 0, 2 = 0.4 and 2 = 1 measured from the sim- 
ulations, together with the linear theory predictions. Linear 
theory predicts that the relative difference of P2 / Po between 
F4 and GR is larger at lower redshifts as expected. How- 
ever, the relative difference measured in simulations on large 
scales is smaller at lower redshifts. This is due to the extra 
damping in the F4 models. This damping becomes stronger 
at lower redshifts and overcomes the linear enhancement of 
f 2 /Po hi the F4 models. 

The strong enhancement in the small scale damping in 
the F4 and F5 models compared to that in GR, which we 
have seen above, could be a clear signal of modified gravity 
that persists at higher redshifts. Of course, here we are only 
talking about the dark matter power spectrum and, as dis- 
cussed in Section 13.21 models for the redshift space power 
spectrum on small scales need to account for nonlinear bias 
effects. As this is most relevant when studying the cluster- 
ing of halos in redshift space we leave this analysis to future 
work. 



4.2 The velocity power spectrum 

The nonlinear evolution of velocity fields on large scales can 
have a significant impact on the redshift space power spec- 
trum. IScoccimarrol (|2004l ) showed that the velocity field is 
more sensitive to tidal gravitational fields compared to the 
density field on large scales. Taking these nonlinear veloc- 
ity effects into account results in an improved model for the 
power spectrum in redshift space i n the quasi-linear regime 
I Jennings. Baugh fc Pascolill2011ah . 

Measuring the velocity power spectrum from sim- 
ulations can be difficult. The method suggested by 
IScoccimarrol (|2004h allows a mass weighted velocity field 
to be constructed but is limited by the fact that it is 
the momentum field which is calculated on a grid and 
so th e velocity field in empty ce lls is artificially set to 
zero ijPueblas fc Scoccimarrol 120091 ). Another limitation of 
this method is that most calculations require the volume 
weighted velocity field instead of the mass weighted field. 



1.4 




z = 1 .0 


1 .2 






1.0 




N 


0.8 




\ 


0.6 







0.01 



0.10 




0.01 



0.10 




0.01 



0.10 



k [h/Mpc] 



Figure 7. The ratio of P2/P0 in the F4 model to that in GR 
at 2 = (bottom panel), 2 = 0.4 (middle panel) and 2 = 1 
(top panel) measured from the simulations. The linear theory 
prediction for this ratio at each redshift is shown as a solid green 
line. 



Using a Delaunay tessellation of a discrete set of points 
allows the desired volume weighted velocity field to be 
constructed accurately on sma ll scales. We use the pub- 
licly available DTFE code (ISchaap fc van de WevgaertJ 



2OOOI; Ivan de Weygaert fc Schaad 2000; 

Cantun fc van de WevgaertJ |201ll ) to construct the ve 



locity divergence field directly. This code constructs the 
Delaunay tessellation from a discrete set of points and 
interpolates the field values onto a user defined grid. For 
the Lbox = 1500 Mpc/h simulation we are able to generate 
the velocity auto, Pee, and cross power spectrum, P$e, on 
a 1024 3 grid. The density field is interpolated onto the grid 
using the cloud-in-cell scheme. The resolution of the mesh 
means that mass assignment effects are negligible on the 
scales of interest here. 
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Figure 8. Left panel: The velocity divergence auto, Peel f 2 (purple) and cross, Pse/f (red) power spectrum at z = measured from 
the F4 model simulation. The lower ratio plot shows Peel Pssf 2 (purple) and Pse/Pssf (red). Right panel: The velocity divergence 
auto, Pee If 2 (blue) and cross, Pse/f (g rev ) power spectrum at z = measured from the ACDM simulation. The lower ratio plot 
shows Peel Pssf 2 (blue) and Pse/Pssf (g rev )- The nonlinear matter power spectrum is plotted as a solid black line for each model. The 
measured P(k) from the L bax = 1500 Mpc/h and L b ox = 1000 Mpc//i simulations are plotted as dashed lines and circles respectively. 



Fig. [8j shows the z = nonlinear velocity and matter 
power spectra measured from the F4 (left panels) and GR 
(right panels) simulations. The errors calculated from the 
scatter amongst six simulations are shown as a hatched re- 
gion for the cross power spectrum, Pse, and as a solid shaded 
region for the auto power spectrum, Pee- We show the ve- 
locity power spectrum from both the L bm = 1500Mpc//i 
(dashed lines) and L box = 1000Mpc//i (circles) simulations. 
The lower panel in each case shows the ratio of Pee /Pssf 2 
and Pse/Pssf for each model. The scales where the lin- 
ear continuity equation breaks down are shown by the de- 
parture of the measured ratios from unity. This occurs on 
slightly larger scales for GR than for the f(R) model. The 
fact that the velocity power spectrum departs from linear 
theory on larger scales than the density field agrees with 
what has been noted previously by IScoccimarr o (2004) and 
Ijennings. Baugh fc Pascolil (|2011al ). We find that Pss and 
Pee differ by up 20% at k < 0.1/i/Mpc. 

When the velocity power spectrum is normalised using 
the linear growth rate as in Fig. [HJ the ratio of the nonlin- 
ear velocity and matter power spectra look very similar in 
both the standard and the modified gravity cosmologies. If 
instead we choose not to normalise 9 using the growth rate, 
we get the curves shown in the top row of Fig.[9]for the ratio 
V Pee /Pss (left panel) and Pse / Pss (right panel) for the F4 
and ACDM cosmologies. Note in linear perturbation theory 
these two ratios equal the linear growth rate, / , which is 
plotted as a solid green line for the f(R) model and as a dot- 
ted black line for ACDM in Fig. 9. There is clearly a large 
difference in the amplitude and shape of these ratios in the 
two models on scales k > 0.03/i/Mpc. It is interesting that 



the ratios agree with the predictions for the linear growth 
rate on scales where the two models can be distinguished 
within the errors, shown by the hatched shaded region for 
f(R) and solid shaded region for ACDM. This is in contrast 
to the multipole moments of the redshift space power spec- 
trum, Fig. U where the Kaiser model predictions using the 
linear growth rate only agree with the measured P(k) on 
extremely large scales where the two cosmologies could not 
be distinguished within the errors. The linear growth rate 
for the F4 model differs from that in ACDM by up to 20% 
for k < 0.1/i/Mpc (see Fig. 1). 

In the middle and bottom rows of Fig [9] we plot similar 
ratios for the F5 and F6 models respectively. For these two 
models the difference in the ratios compared to ACDM is less 
dramatic than for the F4 model. This is to be expected as 
the linear growth rate for these models only differs from that 
in ACDM by at most 6% for F5 and ~ 1% for F6 on scales 
k < O.lft/Mpc. Furthermore, the simulation results start to 
deviate from the linear perturbation prediction earlier than 
it does for the F4 model. This is because the suppression of 
the fifth force itself is a nonlinear effect and the nonlinearity 
gets weaker as \fm\ increases, making the linear perturba- 
tion a better approximation for F4. 

Fig. [10] shows the redshift evolution of the ratio 
\J Pee /Pss = / for the GR (black dotted line for linear per- 
turbation prediction and blue dashed line measured from 
A^-body simulation) and F4 model (green solid and purple 
dashed lines respectively) at z — 0.4 (lower panel) and z = 1 
(upper panel). The absolute difference in this ratio for these 
the two models is even more pronounced at higher redshift, 
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Figure 9. Top left panel: The ratio yj Pgg /Pss at 2 = for the F4 (dashed purple) and ACDM (dashed blue) models. The same ratio 
is shown for F5 and F6 models in the middle and bottom left hand panels respectively. Top right panel: The ratio Psg/Pgs at 2: = for 
the F4 (dashed red) and ACDM (dashed black) models. The same ratio is shown for F5 and F6 in the middle and bottom right hand 
panels respectively. In all panels the variable 6 = V • u. The errors for the f(R) model are shown as hatched shaded regions and as solid 
shaded regions for ACDM. The linear growth rate is plotted as a solid green line and dashed black line for each f(R) model and ACDM 
respectively. 



for the same reason as discussed in the redshift evolution of 
P2/P0 (Fig. H . 

The measured ratio agrees with the linear theory pre- 
dictions for the growth rate on scales k < 0.07/i/Mpc for the 
f(R) model and k < 0.04/i/Mpc for ACDM at z = 1, which 
is again as expected because linear perturbation is a better 
approximation at earlier times. 

Fig. Illl shows the ratio of / in F4 to that in GR at three 
different redshifts. We also plot the linear theory prediction 
for the ratio //(.r)//gr as a green solid line in this figure. 
Linear perturbation theory predicts that the ratio becomes 
larger at lower redshifts, which is shown by a small increase 
in the green line in Fig. 11 at k ~ 1/i/Mpc. On the other 
hand, the ratio //(_r)//gr obtained from the simulations re- 
mains roughly the same at k < 0.2/i/Mpc- at all three red- 



shifts and decreases at increasing redshift on smaller scales 
due to increased damping in the F4 model compared to GR. 
The fractional difference is ~ 12% where the growth rate 
in F4 models peaks and the onset of the increase occurs on 
scales, which is k ~ 0.2,0.15 and 0.09/i/Mpc for 2 = 1,0.4 
and respectively. This is because the damping of the ve- 
locity power spectrum due to nonlinearity becomes larger at 
lower redshifts on small scales, compensating the enhance- 
ment on large scales. The shift of the onset of the peak 
towards larger scales at later times merely reflects the fact 
that small scales are affected earlier. 
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Figure 10. The ratio -J Pee / P55 at z = 0.4 (lower panel) and z = 
1 (upper panel) for the F4 (dashed purple) and ACDM (dashed 
blue) models. The linear growth rate at each rcdshift is plotted as 
a solid green line and dashed black line for the f(R) model and 
ACDM respectively. 



4.3 Extracting the matter and velocity power 
spectra 

In this section we investigate if a model for the 2D redshift 
space power spectrum can be used to extract the density 
and velocity power spectra, as a function of scale, at k < 
O.lft/Mpc. If we were able to measure both of these power 
spectra accurately this would provide us with a measure of 
the growth rate of structure, as seen in Fig. [9] which may be 
scale dependent as is the case for the f{R) gravity model. 
The motivation for restricting ourselves to these large scales, 
k < 0.1/i/Mpc, is that the impact of bias and nonlinear 
damping due to velocit y dispersions is expe cted to be small 
over this range (see e.g. lAngulo et alj|2008h . 

The left and right panels in Fig. [T^] show the 2D power 
spectrum for GR and the F4 model respectively, plotted as 
a function of wavenumber k and fj, at z = 0. The coloured 
shading represent the values of log 10 P(k, /x) measured from 
the simulations. The overplotted red lines represent the 
model of Eq. [25] which uses the velocity and density power 
spectra measured from each simulation. From both of these 
plots it appears that the amplitude predicted by the model 
in Eq. [25] agrees with the measured 2D spectra although it 
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Figure 11. The ratio /f4//gr where / = ^/Pee/Pss at z = 
(bottom panel), z = 0.4 (middle panel) and z = 1.0 (top panel), 
in the F4 model compared to GR (purple dashed lines). The linear 
theory prediction for this ratio at each redshift is shown as a solid 
green line. 



fails to capture the detailed shape of the 2D spectrum over 
the full r ange of k and fi. This r e sult a grees with what was 
found bv lKwan. Lewis "fe Linden (|2012h . 

In order to test the precision with which the model 
in Eq. [25] can reproduce the 2D power spectrum we shall 
fit for both the velocity, Pee, and matter power spectra, 
Pss under the a ssumption that in the q uasi-linear regime 
Pse = V Pee Pss ([Percival fc White! l2009h . We have verified 
that this is true for our simulations to within a few percent 
accuracy for k < 0.1/i/Mpc. We perform this fit over sepa- 
rate k bins of width Ak = 0.01ft/Mpc up to a maximum of 
k — O.lft/Mpc, using the entire range of < fi < 1. 

The results of fitting to the F4 simulation at 2 = and 
z = 1 are shown in the left and right panels respectively in 
Fig. [13] The average power spectra, Pee (lower curves) and 
Pss (upper curves), measured from the six simulations are 
plotted as a black dashed line for ACDM and a solid purple 
line for the f(R) model at each redshift. The red filled circles 
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Figure 12. The two dimensional redshift space power spectrum, P(k,fi), measured from the ACDM simulation (left panel) and the 
F4 simulation (right panel) at z = 0. The coloured contours and black solid lines represent log 10 P(k,fi). The overplottcd red solid lines 
show the predictions of the model of Eq. 1251 where the matter and velocity power spectra used for each cosmology have been measured 
from the simulation. 



show the results of the fit for each power spectra measured 
from the F4 simulation. At z — we show the result of the 
fit for each power spectra measured from the GR simula- 
tion as green squares. At both redshifts it is clear that the 
model in Eq. [25] is able to accurately describe the ampli- 
tude of the 2D P{k,fx). The matter power spectrum, Pss 
is recovered accurately and is distinguishable from ACDM. 
Unfortunately this model is not able to reproduce the ve- 
locity power spectra from the modified gravity model and 
at both redshifts is biased to lower values. We have verified 
that fitting Eq. [25] over a reduced range in (i allows us to 
recover the correct Pee but at the cost of an increase in the 
errors by more than the difference in the two cosmologies. 
These results demonstrate that the redshift space distortion 
model P(k,fi) = P 5S (k) + 2fi 2 y / Pe e Pss + ^ 4 Pee(k), accu- 
rately describes the amplitude of the 2D power spectrum 
and can recover Pss but the angular dependence on \x is in- 
correct and so we canno t extract Pe e - The se results are in 
agreement with work bv lTang et al\ (|201ll ). These authors 
fit for the density velocity cross power spectrum Pse and 
find a similar bias in recovering the velocity power spec- 
trum using this model. We have also fit the model in Eq. 23 
allowing the velocity dispersion damping term to be a free 
parameter. This method recovers the correct Pee within the 
error bars but the constraints on the velocity power spectra 
are too weak to distinguish between GR and the F4 model. 
This measurement of Pee would not be accurate enough to 
allow us to discriminate between the F4 model and ACDM 
at either redshift. 



4.4 Modelling the moments of the redshift space 

P(fc) 

In this section we return to studying the moments of the 
redshift space power spectrum, Po and Pa. As shown in Sec- 
tion [4]3] the model given in Eq. 25 fails to capture the shape 



of the full 2D P(k, fi) so naively we do not expect that we 
can precisely measure these moments on all scales. In this 
section we investigate how well this model works at recov- 
ering the measured moments on large scales after averaging 
over /j,. Previous work has shown that the model given in 
Eq. [2H] provides a good fit to measurements from simula- 
tions on quasi-linear scales k < 0.3/i/Mpc and at high red- 
shifts z ~ 1, without the need to include a damping term 
(|jennings. Baugh fe Pascoli|[2011al ). 



The z = 1 multipole moments, Po (upper curves) and 
P2 (lower curves), measured from the ACDM (empty black 
squares) and F4 (filled purple circles) simulations are shown 
in Fig. 1141 These two power spectra have been separated 
in this plot for clarity. The model given in Eq. [25] using 
the velocity and matter power spectra from the simulations 
are overplotted as a green dot-dashed line and a red dotted 
line for the f(R) and the ACDM cosmologies respectively. 
The inset panel shows the ratio of the measured Po to the 
model for each cosmology, F4 (dot dashed green lines) and 
ACDM (red dotted lines). The model for the monopole mo- 
ment reproduces the measurement for the f(R) model to 
within 10% accuracy at fc < 0.2/i/Mpc. The Po model for 
ACDM is accurate to within 5% at fc < 0.2/i/Mpc. The solid 
black line in the inset panel shows the ratio of the monopole 
moment in the F4 model to ACDM. The model precision 
for ACDM is sufficient to detect the 15% difference in the 
monopole moment which we find between the two cosmolo- 
gies on these large scales. The model in Eq. [26] requires ac- 
curate knowledge of the velocity and matter power spectra 
as input parameters. In this Section we have used Pee and 
Pss measured from the simulations. An alternative to this 
would be to use fitting formula for each of these power spec- 
tra which have sufficient accuracy on these large scales ( see 
e.g. lSmith et aZ.ll2003l ; [Jennings. Baugh fc Pascolill2011ah . 
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Figure 13. Left panel: The nonlinear matter (upper curves) and velocity (lower curves) power spectra measured from the z = F4 
(solid purple) and ACDM (dashed black) simulations. The filled red circles (green squares) show the results from fitting Eq. 25 to the 
2D power spectrum, P(k,fi) from the F4 (GR) simulation. The error bars represent the 1-cr errors on the fit, solid grey and hatched blue 
shaded regions represent the errors on the power spectrum measured from the ACDM and f(R) simulation respectively. Right panel: 
Similar to the left panel but for 2 = 1. 




Figure 14. Upper curves: The monopole moment of the redshift space power spectrum measured from the F4 (filled purple circles) 
and ACDM (empty black squares) simulations at z — 1. Lower curves: The quadrupole moment of the redshift space power spectrum 
measured from the F4 (filled purple circles) and ACDM (empty black squares) simulations at z = 1. The two moments have been offset 
for clarity in this plot. The model in Eg. 1251 for each power spectra moment is shown as a dot-dashed green line for the f(R) model and 
as a red dotted line for ACDM. The hatched (solid) shaded regions represent the errors on the measured power spectra for the f(R) 
(ACDM) simulations. The inset panel shows the ratio of the monopole moment to the model in each cosmology, F4 (dot dashed green 
lines) and ACDM (red dotted lines). The solid black line shows the ratio of the monopole moment in the F4 model to ACDM. 



5 SUMMARY AND CONCLUSIONS 

Modified gravity theories generally predict different cluster- 
ing properties of matter, and as a result both the density 
and the velocity power spectra could be very different from 
the predictions of general relativity. The f(R) gravity model 



has been a leading example of this in recent years. Here, the 
enhancement to the standard gravity depends sensitively on 
the local matter density through the so-called chameleon 
mechanism. In high matter density and high curvature re- 
gions (fa <C |$| where <£> is the Newtonian potential), the 
enhancement is strongly suppressed and standard gravity is 
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recovered; on the other hand, in low matter density and low 
curvature regions the enhancement factor can be as large 
as 4/3. Depending on the value of /ho and the local en- 
vironment, the transition scale, or Compton length, of the 
scalaron ranges from less than one to more than a few mega- 
parsecs, which can potentially leave detectable features in 
the distribution of matter and clustering patterns of the 
large-scale structure. 

Galaxy surveys measure the distribution of matter in 
redshift space where the true position of a galaxy appears 
distorted along the line of sight due to peculiar velocities. 
The goal of many current and furture galaxy redshift surveys 
is to constrain deviations from GR and so it is important 
to understand how observables are affected by redshift space 
distortions. Theoretical studies of this require both high res- 
olution and large volume numerical simulations, which pre- 
viously have not been performed for modified gravity mod- 
els. In this paper we use simulations in large volume boxes to 
carry out the first study of the clustering of the dark matter 
in redshift space in a f(R) modifie d gravity. The simulation 
code developed bv lLi et al\ l|2012h allows us to model large 
volumes with good resolution. 

The simulations we use in this analysis have two differ- 
ent resolutions, with 1024 3 dark matter particles in compu- 
tational boxes of Lbox = 1.5/i _1 Gpc and 1.0/i _1 Gpc on a 
side. We have checked that the simulations agree with one 
another down to k ~ 0.3/i/Mpc. We have compared the mat- 
ter P(k) measured in real space against previous simulations 
and found good agreement. 

We have measured the redshift space power spectrum in 
the GR and the f(R) cosmologies at redshifts 2 = 0, 0.4 and 
1. We find an enhanced boost in the power on large scales 
and a substantial increase in the damping on small scales in 
the f(R) cosmology compared to GR at all redshifts. The de- 
viations are largest for the f(R) model with parameter value 
\fno\ = 10 -4 (F4) and are reduced as \fno\ is decreased. 
The large scale enhancement of the power is a result of the 
fifth force in the modified gravity theory, which strengthens 
the matter clustering on large scales. On small scales where 
the local curvature is not too high, the fifth force makes 
the particles move faster, increasing the velocity dispersion 
and causing a stronger damping of the power compared to 
GR. These effects can be seen at various redshifts, and are 
even slightly stronger at earlier times. However, for some 
model parameters such as |/ro| = 10~ 6 , the fifth force is 
strongly suppressed by the chameleon mechanism, and its 
effect is too weak to be distinguished from standard gravity. 
We find that the relative difference in the moments of the 
power spectrum, P2/P0, between the F4 f(R) model and GR 
ranges from 20% at 2 = 1 to 40% at 2 = at k = 0.2/i/Mpc 
due to the enhanced nonlinear damping on small scales. 

We measure the velocity divergence power spectrum us- 
ing the DTFE method in both f(R) and GR cosmologies 
and find a large difference between Peg in the f(R) model 
compared to GR. This difference is much larger than the 
difference between the non linear matter P(k) in the two 
models suggesting that the velocity power spectra is a far 
more sensitive probe of modified gravity. We find a large 
deviation in the ratios y/ Pee / Pss and Pse /Pss between the 
two models at 0.03 < fe(/i/Mpc)< 0.5 at z = 0. In linear 
theory these ratios equal the growth rate of structure, /, 
on large scales when the velocity divergence is normalized 



as 9 — V ■ v/aH. Our results show that the measured ra- 
tios agree with the linear growth rate for each cosmology, 
which is scale dependent in the case of modified gravity, for 
k < 0.06/i/Mpc at z = 0. We find that the relative deviation 
of the measured nonlinear ratio, \J Pee / Pss , from the linear 
prediction between the f(R) model and GR decreases with 
increasing redshift. 

Using a simple quasi linear model for the 2D redshift 
space power spectrum, Eq. Q25p . which includes nonlinear 
velocity terms but no small scale damping parameter, we 
attempt to extract the matter and velocity power spectra. 
On scales k < 0.1/i/Mpc we can recover the nonlinear matter 
power spectrum to within a few percent for both the f(R) 
and the GR cosmology. The model fails to describe the shape 
of the 2D power spectrum and we are unable to reconstruct 
the velocity P(k) accurately. The fact that this model re- 
covers the non linear matter power spectrum so precisely 
indicates that this method can be used to contrain modified 
gravity models. Our simulation results show that improved 
theoretical models are required in order to measure the ve- 
locity power spectrum where there is a large difference in 
the predicted signal between these two cosmologies. 

We show that the same model works very well at fitting 
the first two multipole moments of the redshift space power 
spectrum on large scales, especially at high redshifts. We are 
able to match the monopole moment to within 3% for GR 
and 10% for the f(R) cosmology on scales k < 0.2/i/Mpc at 
2 = 1. This difference is smaller than the 15% difference in 
the F4 and GR models on the same scales. 

In this study we have addressed two separate questions. 
The first considers how well a simple model for redshift space 
distortions, Eq. 25, works at recovering the nonlinear mat- 
ter and velocity P(k) and how well it describes the mul- 
tipole moments of the redshift space power spectrum. The 
matter P(k) and the monopole moment, Po are recovered 
accurately and as a result we do not think that the failure 
of the model to decribe the full 2D shape of P(fc,/i) is a 
serious pitfall for future redshift space analyses. Our results 
point to necessary adjustments needed in the model, such as 
including non-linear damping terms for example. Whether 
the velocity power spectrum can be extracted using more 
complicated modelling of the 2D power spectrum, and how 
well it can be measured, is left for future analysis. The sec- 
ond, two-part, question is which observable, P(k,^), Pq, P2 
etc. shows the largest difference between an f(R) and a GR 
cosmology and is the simple redshift space distortion model 
accurate enough to allow us to measure these differences. We 
observe the largest difference between these two cosmolo- 
gies in the measured non linear velocity power spectrum on 
scales k > 0.03h/Mpc and in the ratio of the multipole mo- 
ments P2/P0 on scales k > Q.2h/Mpc. Our results show that 
an improved model for redshift space distortions is needed 
in order to extract the velocity power spectrum to a suffi- 
cient accuracy to distinguish a f(R) from a GR cosmology. 
The differences between the ratios of the multipole moments, 
P2/P0, in GR and f(R) are mainly on small scales. We shall 
address whether or not this difference is present in the red- 
shift space clustering of halos and the effects of bias in future 
work. 

To conclude, we find that redshift space distortions in 
modified gravity models have an impact on the clustering of 
dark matter on large and small scales to a level which may be 
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distinguished from general relativity. Current redshift space 
distortion models, which are valid on quasi linear scales, are 
accurate enough to extract the non linear matter P(k) in real 
space from the measured 2D redshift space power spectrum 
on large scales, allowing us to constrain f(R) modified grav- 
ity. The large difference between the predicted velocity P(k) 
in the f(R) and GR cosmology make this a very promising 
observable with which to test GR provided that this can 
be accurately extracted from the redshift space power spec- 
trum. For certain f(R) parameter values, e.g. |/ro| = 1CP 6 , 
the impact of modified gravity on both the matter and ve- 
locity fields is not significant and any deviations from GR 
are restricted to small scales. In a follow up paper we will 
examine the redshift space distortions in the clustering of 
halos in these modified gravity models on nonlinear scales. 
We will test several nonlinear models for redshift space dis- 
tortions to predict how well this observable can constrain 
f(R) modified gravity in future surveys. 
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